Dependencies

library(ggridges) # ggarrange
library(ggpubr)
library(ggdist) 
library(ggtext)
library(ggforce)

library(tidyverse) # data manipulation
library(modelr) # data manipulation

library(rstan) #detectCores
library(brms) # core pack
library(parallel) #detectCores
library(tidybayes) 

# a flag to let us load models from the disk (F) or re-run the models (T)
# it took about 20 ~ 60 minutes to run the models
REBUILD_FLAG <<- FALSE

Load data

metas <- read.csv("../data/metas.csv", stringsAsFactors = F, header = T)
trials_all <- read.csv("../data/trials.csv", stringsAsFactors = F, header = T)
seper <- " o "

trials_all <- trials_all %>%
          filter(TrialStatus == "main") %>% #drop practice trial
          mutate(Error = TrialAnswer - TrialIntAccuracy,
                 ErrorMagnitude = abs(Error),
                 ResponseTime = (TrialVisualizationEnd - TrialVisualizationStart) / 1000,
                 LogResponseTime = log(ResponseTime)) %>% # log transform
          filter(!(ParticipantID == "P7" 
                   & SessionIndex == 0 
                   & TrialStatus == "main"
                   & TrialIndex == 16)) %>%  # filter mistakes
          filter(!(ParticipantID == "P9" 
                   & SessionIndex == 2 
                   & TrialStatus == "main" 
                   & TrialIndex == 17)) %>% 
          filter(!(ParticipantID == "P11" 
                   & SessionIndex == 0 
                   & TrialStatus == "main"  
                   & TrialIndex == 12)) %>% 
          filter(!(ParticipantID == "P20" 
                   & SessionIndex == 3 
                   & TrialStatus == "main" 
                   & TrialIndex == 5)) %>% 
          filter(!(ParticipantID == "P24" 
                   & SessionIndex == 0 
                   & TrialStatus == "main"  
                   & TrialIndex == 26)) %>% 
          filter(ErrorMagnitude < 51) %>% 
          mutate(  
                 # manual recode levels using orthogonal contrast coding  
                 # the distance within a variable is always 1
                 AllFactors = paste(FactorStereo, FactorMotion, FactorPerspective, FactorShading, Dim, sep = seper),
                 Device = recode(FactorStereo, yes = 0.5, no = -0.5),
                 Motion = recode(FactorMotion, yes = 0.5, no = -0.5),
                 Projection = recode(FactorPerspective, yes = 0.5, no = -0.5),
                 # this will ignore shading-1, the commonly-used phong lighting model
                 Shading = recode(FactorShading, `shading-1` = 0, `shading-2` = 0.5, `flat`= -0.5, `no` = -0.5),
                 DataModel = recode(TrialDataset, "cifar10" = 0.5, "babi-q1" = -0.5),
                 Dimensionality = recode(Dim, `2` = -0.5, `3` = 0.5)
                 ) 


trials <- trials_all %>%
  select(ParticipantID, TrialIntAccuracy, ErrorMagnitude, ResponseTime, Device, Motion, Projection, Shading, Dimensionality, DataModel)


vtrials <- trials_all %>%
  select(ParticipantID, ErrorMagnitude, ResponseTime, FactorStereo, FactorMotion, FactorPerspective, FactorShading, Dim, TrialDataset, AllFactors)

Raw data

We first plot all raw data to get a sense of how data distribute. Both error magnitude and response time datum seem not to follow a normal distribution, and they are both non-negative. A skewed normal, lognormal, and gamma distribution could be reasonable. Note that error magnitude contains 0, and response time has a upper bound 60s (e.g., 60.1s, 60.2s due to running time). We will use gamma distribution for both measures. To formally compare which distribution better fits data, we could use the model_diagnostics function defined below.

g_raw_error <- vtrials %>%
         ggplot(mapping = aes(x = ErrorMagnitude)) +
         facet_wrap(AllFactors ~ .,ncol = 3, scales = "free") +
         coord_cartesian(xlim = c(0, 50), ylim = c(0, 30)) + 
         geom_histogram(binwidth = 2, fill = "lightblue", color = "white") +
         ggtitle("raw data - error magnitude (pt)") + 
         xlab ("error magnitude ") + ylab("count")

ggsave(g_raw_error, filename = "pic/g_raw_error.pdf", width = 10, height = 10, units = "in")
g_raw_error

g_raw_time <- vtrials %>%
         ggplot(mapping = aes(x = ResponseTime)) +
         facet_wrap(AllFactors ~ .,ncol = 3, scales = "free") +
         coord_cartesian(xlim = c(0, 62), ylim = c(0, 30)) + 
         geom_histogram(binwidth = 2, fill = "lightgray", color = "white") +
         ggtitle("raw data - response time  (s)") + 
         xlab ("response time ") + ylab("count")

ggsave(g_raw_time, filename = "pic/g_raw_time.pdf", width = 10, height = 10, units = "in")

g_raw_time

# Models

Helper function

First we define a helper function to wrap up model diagnostics.

model_diagnostics <- function(m, name = "model diagnostics"){
  print(name)
  
  print(summary(m))

  # plot(m) # run this line to check convergence when first fitting the model
  
  plot(pp_check(m, type = "dens_overlay", nsamples = 100))
  plot(pp_check(m, type = "stat_2d", nsamples = 2000))
  plot(pp_check(m, type = "error_binned", nsamples = 100))
  
  print(loo(m))
  print(waic(m))
}

Model specification

For error magnitude, we have 0. Therefore, we use a hurdle_gamma distribution to account for this. Based on the experimental design, Dimensionality only interacts with Device and DataModel. The other five variables could interact with each other. We are also concerned about the difference between Dimensionality; we therefore model it as a group-level effect. The two data models have different ranges of accuracy; we also specify it as a group-effect. Last, each trial has different true accuracies, and we specify it as another random intercept to separate the effects. Because of the skewness (we also don’t have evidence of error magnitude being normally distributed), we use a gamma likelihood instead of normal/student_t. Because of the zero values (no error), we use a hurdle gamma likelihood. Other choices like skewed normal are also reasonable, but gamma results in posterior prediction best capturing the data and of a smaller waic. Since we care about mean of error magnitude and the data also does not suggest a big variance in standard deviation across conditions, we did not use a submodel for standard deviation.

if(REBUILD_FLAG){

  m_error_hurdle_gamma <- bf(ErrorMagnitude ~  
                           0 + Intercept + 
                           Device*Motion*Shading*Projection*DataModel +  
                           Device*Dimensionality*DataModel + 
                           (1 + Dimensionality*DataModel | ParticipantID) + 
                           (1 | TrialIntAccuracy), 
                           family = hurdle_gamma(link = "log")
                        ) 
  
  get_prior(m_error_hurdle_gamma, data = trials)

    
  m_brms_error_hurdle_gamma <<- brm(m_error_hurdle_gamma, 
                               data = trials, 
                               family = hurdle_gamma(link = "log"), 
                               prior = c( # priors" 2 SDs roughly catch 95% data
                                  prior(normal(0, 2.5), class = "b"),
                                  prior(normal(0, 25),  class = "b", coef = "Intercept"), 
                                  prior(gamma(2, 10), class = "shape"),
                                  prior(student_t(3, 0, 2.5), class = "sd")
                                ), 
                                backend = 'cmdstanr',
                                chains = N_CHAINS, iter = N_ITER, warmup = N_WARMUP, thin = N_THIN, cores = N_CORES,
                                control = list(adapt_delta = DELTA, max_treedepth = TREE_DEPTH))
  
  saveRDS(file = "rds/m_brms_error_hurdle_gamma.rds", m_brms_error_hurdle_gamma)
  
}else{
  m_brms_error_hurdle_gamma <<- readRDS("rds/m_brms_error_hurdle_gamma.rds")
}

model_diagnostics(m_brms_error_hurdle_gamma)
## [1] "model diagnostics"
##  Family: hurdle_gamma 
##   Links: mu = log; shape = identity; hu = identity 
## Formula: ErrorMagnitude ~ 0 + Intercept + Device * Motion * Shading * Projection * DataModel + Device * Dimensionality * DataModel + (1 + Dimensionality * DataModel | ParticipantID) + (1 | TrialIntAccuracy) 
##    Data: trials (Number of observations: 3320) 
##   Draws: 2 chains, each with iter = 2000; warmup = 0; thin = 2;
##          total post-warmup draws = 2000
## 
## Group-Level Effects: 
## ~ParticipantID (Number of levels: 32) 
##                                              Estimate Est.Error l-95% CI
## sd(Intercept)                                    0.21      0.04     0.14
## sd(Dimensionality)                               0.14      0.08     0.01
## sd(DataModel)                                    0.38      0.07     0.25
## sd(Dimensionality:DataModel)                     0.20      0.15     0.01
## cor(Intercept,Dimensionality)                    0.15      0.37    -0.58
## cor(Intercept,DataModel)                        -0.16      0.23    -0.59
## cor(Dimensionality,DataModel)                   -0.35      0.34    -0.88
## cor(Intercept,Dimensionality:DataModel)         -0.08      0.41    -0.81
## cor(Dimensionality,Dimensionality:DataModel)    -0.13      0.43    -0.86
## cor(DataModel,Dimensionality:DataModel)          0.10      0.41    -0.72
##                                              u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                                    0.30 1.00     2461     3523
## sd(Dimensionality)                               0.30 1.00     1341     2692
## sd(DataModel)                                    0.54 1.00     2566     3284
## sd(Dimensionality:DataModel)                     0.54 1.00     1830     2654
## cor(Intercept,Dimensionality)                    0.82 1.00     3101     3255
## cor(Intercept,DataModel)                         0.32 1.00     1920     2239
## cor(Dimensionality,DataModel)                    0.40 1.00      661     1023
## cor(Intercept,Dimensionality:DataModel)          0.72 1.00     3422     3497
## cor(Dimensionality,Dimensionality:DataModel)     0.73 1.00     2625     3253
## cor(DataModel,Dimensionality:DataModel)          0.83 1.00     3607     3287
## 
## ~TrialIntAccuracy (Number of levels: 45) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.31      0.04     0.24     0.40 1.00     1935     2886
## 
## Population-Level Effects: 
##                                            Estimate Est.Error l-95% CI u-95% CI
## Intercept                                      2.02      0.07     1.88     2.16
## Device                                        -0.01      0.07    -0.15     0.12
## Motion                                        -0.14      0.04    -0.21    -0.06
## Shading                                       -0.03      0.03    -0.10     0.04
## Projection                                    -0.01      0.03    -0.06     0.05
## DataModel                                      0.06      0.10    -0.14     0.25
## Dimensionality                                 0.10      0.07    -0.04     0.25
## Device:Motion                                  0.01      0.07    -0.14     0.16
## Device:Shading                                -0.01      0.07    -0.15     0.12
## Motion:Shading                                 0.08      0.07    -0.05     0.22
## Device:Projection                             -0.03      0.06    -0.15     0.08
## Motion:Projection                             -0.04      0.06    -0.15     0.07
## Shading:Projection                             0.06      0.07    -0.08     0.19
## Device:DataModel                               0.12      0.14    -0.15     0.38
## Motion:DataModel                              -0.09      0.07    -0.24     0.05
## Shading:DataModel                             -0.05      0.07    -0.19     0.09
## Projection:DataModel                           0.08      0.06    -0.03     0.20
## Device:Dimensionality                          0.13      0.14    -0.13     0.41
## DataModel:Dimensionality                       0.13      0.14    -0.16     0.40
## Device:Motion:Shading                          0.07      0.14    -0.21     0.34
## Device:Motion:Projection                       0.07      0.11    -0.14     0.29
## Device:Shading:Projection                      0.22      0.14    -0.04     0.48
## Motion:Shading:Projection                      0.04      0.14    -0.23     0.33
## Device:Motion:DataModel                        0.06      0.15    -0.23     0.35
## Device:Shading:DataModel                       0.02      0.14    -0.25     0.29
## Motion:Shading:DataModel                       0.40      0.14     0.12     0.66
## Device:Projection:DataModel                    0.14      0.11    -0.07     0.36
## Motion:Projection:DataModel                   -0.01      0.11    -0.23     0.21
## Shading:Projection:DataModel                  -0.12      0.14    -0.40     0.15
## Device:DataModel:Dimensionality               -0.12      0.27    -0.65     0.40
## Device:Motion:Shading:Projection               0.00      0.28    -0.53     0.55
## Device:Motion:Shading:DataModel                0.24      0.26    -0.27     0.76
## Device:Motion:Projection:DataModel            -0.20      0.23    -0.65     0.25
## Device:Shading:Projection:DataModel           -0.01      0.27    -0.54     0.53
## Motion:Shading:Projection:DataModel            0.45      0.27    -0.10     0.98
## Device:Motion:Shading:Projection:DataModel    -0.10      0.55    -1.19     0.95
##                                            Rhat Bulk_ESS Tail_ESS
## Intercept                                  1.00     1954     3030
## Device                                     1.00     3520     3500
## Motion                                     1.00     3542     3488
## Shading                                    1.00     3429     3319
## Projection                                 1.00     3491     3455
## DataModel                                  1.00     3159     3204
## Dimensionality                             1.00     3189     3462
## Device:Motion                              1.00     3465     3294
## Device:Shading                             1.00     3442     3316
## Motion:Shading                             1.00     3552     2685
## Device:Projection                          1.00     3244     3246
## Motion:Projection                          1.00     3258     3183
## Shading:Projection                         1.00     3360     3312
## Device:DataModel                           1.00     3770     3464
## Motion:DataModel                           1.00     3505     3692
## Shading:DataModel                          1.00     3368     3439
## Projection:DataModel                       1.00     3441     2814
## Device:Dimensionality                      1.00     3593     3604
## DataModel:Dimensionality                   1.00     3131     3450
## Device:Motion:Shading                      1.00     3595     3334
## Device:Motion:Projection                   1.00     3551     3228
## Device:Shading:Projection                  1.00     3297     3229
## Motion:Shading:Projection                  1.00     3468     3099
## Device:Motion:DataModel                    1.00     3667     3435
## Device:Shading:DataModel                   1.00     3470     3180
## Motion:Shading:DataModel                   1.00     2700     3165
## Device:Projection:DataModel                1.00     3083     2932
## Motion:Projection:DataModel                1.00     3354     3284
## Shading:Projection:DataModel               1.00     3367     3183
## Device:DataModel:Dimensionality            1.00     3883     3283
## Device:Motion:Shading:Projection           1.00     3362     3072
## Device:Motion:Shading:DataModel            1.00     3562     3169
## Device:Motion:Projection:DataModel         1.00     3047     3315
## Device:Shading:Projection:DataModel        1.00     3412     3486
## Motion:Shading:Projection:DataModel        1.00     3524     3156
## Device:Motion:Shading:Projection:DataModel 1.00     3376     2850
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     1.80      0.04     1.72     1.89 1.00     3234     2612
## hu        0.07      0.00     0.06     0.08 1.00     3309     3104
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

## 
## Computed from 4000 by 3320 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -9723.0 49.6
## p_loo       129.7  4.6
## looic     19446.0 99.2
## ------
## Monte Carlo SE of elpd_loo is 0.3.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     3316  99.9%   530       
##  (0.5, 0.7]   (ok)          4   0.1%   351       
##    (0.7, 1]   (bad)         0   0.0%   <NA>      
##    (1, Inf)   (very bad)    0   0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
## 
## Computed from 4000 by 3320 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic  -9721.4 49.5
## p_waic       128.2  4.4
## waic       19442.9 99.0
## 
## 32 (1.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

For response time, we have a similar motivation with error magnitude with a few exceptions. First response time is positive. Therefore a gamma distribution is good enough. Second, we know the upper bound of response time is about 60 is cropped. We need to tell the model about this (left censored). Also, data model and dimensionality probably does not affect response time very much; however, different devices and motion levels should be treated differently. So we use Device*Motion as the group-level effects to replace data properties as group-level effects. We note that the distributions of response time could look very different; a sub-model for shape is more appropriate to capture the difference. We added Device and Motion as fixed effects (adding random effects would let the model not converge.)

if(REBUILD_FLAG == TRUE){
  
  trials <- trials %>%
      mutate(cen = ifelse(ResponseTime > 60, "left", "none"))

  m_time_gamma <- bf(
                    ResponseTime | cens(cen) ~ 0 + Intercept +
                      Device*Motion*Shading*Projection*DataModel +  
                      Device*Dimensionality*DataModel + 
                      (1 + Device*Motion | ParticipantID),
                    shape ~ 0 + Intercept + 
                      Device*Motion,
                      family = Gamma(link="log")
                     ) 
    
  priors_time_gamma <- c(
                prior(normal(0, 2.5), class = "b"),
                prior(normal(0, 30), class = "b", coef = "Intercept"), 
                #prior(gamma(1, 20), class = "shape"),
                prior(student_t(3, 0, 2.5), class = "sd")
                ) 
    
  m_brms_time_gamma <<- brm(m_time_gamma, 
                              data = trials, family = Gamma(link="log"),  prior = priors_time_gamma, 
                              backend = 'cmdstanr',
                              chains = N_CHAINS, iter = N_ITER, warmup = N_WARMUP, thin = N_THIN, cores = N_CORES,
                              control = list(adapt_delta = DELTA, max_treedepth = TREE_DEPTH))
  
  saveRDS(file = "rds/m_brms_time_gamma.rds", m_brms_time_gamma)
  
}else{
  
  m_brms_time_gamma <<- readRDS("rds/m_brms_time_gamma.rds")
}

model_diagnostics(m_brms_time_gamma)
## [1] "model diagnostics"
##  Family: gamma 
##   Links: mu = log; shape = log 
## Formula: ResponseTime | cens(cen) ~ 0 + Intercept + Device * Motion * Shading * Projection * DataModel + Device * Dimensionality * DataModel + (1 + Device * Motion | ParticipantID) 
##          shape ~ 0 + Intercept + Device * Motion
##    Data: trials (Number of observations: 3320) 
##   Draws: 2 chains, each with iter = 2000; warmup = 0; thin = 2;
##          total post-warmup draws = 2000
## 
## Group-Level Effects: 
## ~ParticipantID (Number of levels: 32) 
##                              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                    0.64      0.08     0.50     0.83 1.00     1523
## sd(Device)                       0.25      0.04     0.18     0.34 1.00     2253
## sd(Motion)                       0.29      0.05     0.21     0.39 1.00     2545
## sd(Device:Motion)                0.75      0.11     0.56     1.00 1.00     2499
## cor(Intercept,Device)           -0.11      0.19    -0.46     0.26 1.00     3245
## cor(Intercept,Motion)           -0.48      0.15    -0.73    -0.14 1.00     2975
## cor(Device,Motion)              -0.10      0.20    -0.48     0.29 1.00     2583
## cor(Intercept,Device:Motion)    -0.14      0.18    -0.47     0.22 1.00     2874
## cor(Device,Device:Motion)        0.05      0.20    -0.36     0.42 1.00     2015
## cor(Motion,Device:Motion)        0.08      0.19    -0.30     0.44 1.00     1829
##                              Tail_ESS
## sd(Intercept)                    2186
## sd(Device)                       2867
## sd(Motion)                       3241
## sd(Device:Motion)                2510
## cor(Intercept,Device)            3523
## cor(Intercept,Motion)            3479
## cor(Device,Motion)               3367
## cor(Intercept,Device:Motion)     3201
## cor(Device,Device:Motion)        2950
## cor(Motion,Device:Motion)        2740
## 
## Population-Level Effects: 
##                                            Estimate Est.Error l-95% CI u-95% CI
## Intercept                                      2.27      0.11     2.05     2.50
## Device                                         0.37      0.06     0.25     0.50
## Motion                                         0.59      0.05     0.49     0.70
## Shading                                        0.01      0.02    -0.04     0.05
## Projection                                     0.04      0.02     0.00     0.08
## DataModel                                      0.14      0.07    -0.00     0.28
## Dimensionality                                 0.12      0.05     0.02     0.21
## Device:Motion                                  0.37      0.14     0.10     0.65
## Device:Shading                                 0.09      0.05    -0.00     0.18
## Motion:Shading                                 0.14      0.05     0.05     0.24
## Device:Projection                              0.03      0.04    -0.04     0.11
## Motion:Projection                              0.05      0.04    -0.02     0.13
## Shading:Projection                            -0.07      0.05    -0.17     0.02
## Device:DataModel                               0.03      0.14    -0.26     0.31
## Motion:DataModel                              -0.02      0.13    -0.29     0.23
## Shading:DataModel                              0.07      0.05    -0.02     0.17
## Projection:DataModel                           0.02      0.04    -0.06     0.09
## Device:Dimensionality                         -0.06      0.10    -0.26     0.13
## DataModel:Dimensionality                       0.05      0.10    -0.13     0.25
## Device:Motion:Shading                         -0.06      0.09    -0.24     0.12
## Device:Motion:Projection                      -0.05      0.08    -0.20     0.10
## Device:Shading:Projection                      0.07      0.09    -0.11     0.26
## Motion:Shading:Projection                      0.10      0.10    -0.09     0.29
## Device:Motion:DataModel                        0.39      0.26    -0.11     0.90
## Device:Shading:DataModel                      -0.02      0.10    -0.20     0.17
## Motion:Shading:DataModel                       0.05      0.09    -0.13     0.24
## Device:Projection:DataModel                   -0.03      0.08    -0.18     0.13
## Motion:Projection:DataModel                   -0.06      0.08    -0.21     0.10
## Shading:Projection:DataModel                   0.01      0.09    -0.17     0.19
## Device:DataModel:Dimensionality               -0.06      0.19    -0.43     0.31
## Device:Motion:Shading:Projection               0.24      0.18    -0.13     0.59
## Device:Motion:Shading:DataModel               -0.12      0.19    -0.48     0.24
## Device:Motion:Projection:DataModel             0.08      0.15    -0.22     0.38
## Device:Shading:Projection:DataModel           -0.05      0.19    -0.42     0.32
## Motion:Shading:Projection:DataModel            0.27      0.19    -0.10     0.63
## Device:Motion:Shading:Projection:DataModel     0.20      0.37    -0.51     0.93
## shape_Intercept                                1.30      0.02     1.26     1.35
## shape_Device                                   0.22      0.05     0.13     0.31
## shape_Motion                                   0.32      0.05     0.22     0.41
## shape_Device:Motion                            0.20      0.10     0.01     0.40
##                                            Rhat Bulk_ESS Tail_ESS
## Intercept                                  1.00     1361     2099
## Device                                     1.00     3072     3488
## Motion                                     1.00     2404     3262
## Shading                                    1.00     3621     3624
## Projection                                 1.00     3976     3552
## DataModel                                  1.00     2930     3508
## Dimensionality                             1.00     3571     3331
## Device:Motion                              1.00     2564     3077
## Device:Shading                             1.00     3496     3657
## Motion:Shading                             1.00     3564     3428
## Device:Projection                          1.00     3704     3397
## Motion:Projection                          1.00     3707     3491
## Shading:Projection                         1.00     3915     3561
## Device:DataModel                           1.00     3009     3055
## Motion:DataModel                           1.00     2601     2912
## Shading:DataModel                          1.00     3854     3575
## Projection:DataModel                       1.00     3657     3362
## Device:Dimensionality                      1.00     3402     3159
## DataModel:Dimensionality                   1.00     3372     3664
## Device:Motion:Shading                      1.00     3559     3334
## Device:Motion:Projection                   1.00     3971     3807
## Device:Shading:Projection                  1.00     3918     3685
## Motion:Shading:Projection                  1.00     3751     3604
## Device:Motion:DataModel                    1.00     2701     3242
## Device:Shading:DataModel                   1.00     3919     3614
## Motion:Shading:DataModel                   1.00     3616     3661
## Device:Projection:DataModel                1.00     3600     3452
## Motion:Projection:DataModel                1.00     3756     3346
## Shading:Projection:DataModel               1.00     3820     3572
## Device:DataModel:Dimensionality            1.00     3543     3879
## Device:Motion:Shading:Projection           1.00     3723     3436
## Device:Motion:Shading:DataModel            1.00     3533     3542
## Device:Motion:Projection:DataModel         1.00     3781     3634
## Device:Shading:Projection:DataModel        1.00     3898     3767
## Motion:Shading:Projection:DataModel        1.00     3679     3492
## Device:Motion:Shading:Projection:DataModel 1.00     3725     3657
## shape_Intercept                            1.00     3592     3524
## shape_Device                               1.00     3810     3379
## shape_Motion                               1.00     3590     3237
## shape_Device:Motion                        1.00     3736     3499
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

## 
## Computed from 4000 by 3320 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo  -9796.9  74.5
## p_loo       179.8  10.5
## looic     19593.9 149.0
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     3301  99.4%   555       
##  (0.5, 0.7]   (ok)         15   0.5%   155       
##    (0.7, 1]   (bad)         4   0.1%   37        
##    (1, Inf)   (very bad)    0   0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.
## 
## Computed from 4000 by 3320 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic  -9795.1  74.3
## p_waic       178.0  10.2
## waic       19590.1 148.7
## 
## 68 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.
summary(m_brms_time_gamma)
##  Family: gamma 
##   Links: mu = log; shape = log 
## Formula: ResponseTime | cens(cen) ~ 0 + Intercept + Device * Motion * Shading * Projection * DataModel + Device * Dimensionality * DataModel + (1 + Device * Motion | ParticipantID) 
##          shape ~ 0 + Intercept + Device * Motion
##    Data: trials (Number of observations: 3320) 
##   Draws: 2 chains, each with iter = 2000; warmup = 0; thin = 2;
##          total post-warmup draws = 2000
## 
## Group-Level Effects: 
## ~ParticipantID (Number of levels: 32) 
##                              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                    0.64      0.08     0.50     0.83 1.00     1523
## sd(Device)                       0.25      0.04     0.18     0.34 1.00     2253
## sd(Motion)                       0.29      0.05     0.21     0.39 1.00     2545
## sd(Device:Motion)                0.75      0.11     0.56     1.00 1.00     2499
## cor(Intercept,Device)           -0.11      0.19    -0.46     0.26 1.00     3245
## cor(Intercept,Motion)           -0.48      0.15    -0.73    -0.14 1.00     2975
## cor(Device,Motion)              -0.10      0.20    -0.48     0.29 1.00     2583
## cor(Intercept,Device:Motion)    -0.14      0.18    -0.47     0.22 1.00     2874
## cor(Device,Device:Motion)        0.05      0.20    -0.36     0.42 1.00     2015
## cor(Motion,Device:Motion)        0.08      0.19    -0.30     0.44 1.00     1829
##                              Tail_ESS
## sd(Intercept)                    2186
## sd(Device)                       2867
## sd(Motion)                       3241
## sd(Device:Motion)                2510
## cor(Intercept,Device)            3523
## cor(Intercept,Motion)            3479
## cor(Device,Motion)               3367
## cor(Intercept,Device:Motion)     3201
## cor(Device,Device:Motion)        2950
## cor(Motion,Device:Motion)        2740
## 
## Population-Level Effects: 
##                                            Estimate Est.Error l-95% CI u-95% CI
## Intercept                                      2.27      0.11     2.05     2.50
## Device                                         0.37      0.06     0.25     0.50
## Motion                                         0.59      0.05     0.49     0.70
## Shading                                        0.01      0.02    -0.04     0.05
## Projection                                     0.04      0.02     0.00     0.08
## DataModel                                      0.14      0.07    -0.00     0.28
## Dimensionality                                 0.12      0.05     0.02     0.21
## Device:Motion                                  0.37      0.14     0.10     0.65
## Device:Shading                                 0.09      0.05    -0.00     0.18
## Motion:Shading                                 0.14      0.05     0.05     0.24
## Device:Projection                              0.03      0.04    -0.04     0.11
## Motion:Projection                              0.05      0.04    -0.02     0.13
## Shading:Projection                            -0.07      0.05    -0.17     0.02
## Device:DataModel                               0.03      0.14    -0.26     0.31
## Motion:DataModel                              -0.02      0.13    -0.29     0.23
## Shading:DataModel                              0.07      0.05    -0.02     0.17
## Projection:DataModel                           0.02      0.04    -0.06     0.09
## Device:Dimensionality                         -0.06      0.10    -0.26     0.13
## DataModel:Dimensionality                       0.05      0.10    -0.13     0.25
## Device:Motion:Shading                         -0.06      0.09    -0.24     0.12
## Device:Motion:Projection                      -0.05      0.08    -0.20     0.10
## Device:Shading:Projection                      0.07      0.09    -0.11     0.26
## Motion:Shading:Projection                      0.10      0.10    -0.09     0.29
## Device:Motion:DataModel                        0.39      0.26    -0.11     0.90
## Device:Shading:DataModel                      -0.02      0.10    -0.20     0.17
## Motion:Shading:DataModel                       0.05      0.09    -0.13     0.24
## Device:Projection:DataModel                   -0.03      0.08    -0.18     0.13
## Motion:Projection:DataModel                   -0.06      0.08    -0.21     0.10
## Shading:Projection:DataModel                   0.01      0.09    -0.17     0.19
## Device:DataModel:Dimensionality               -0.06      0.19    -0.43     0.31
## Device:Motion:Shading:Projection               0.24      0.18    -0.13     0.59
## Device:Motion:Shading:DataModel               -0.12      0.19    -0.48     0.24
## Device:Motion:Projection:DataModel             0.08      0.15    -0.22     0.38
## Device:Shading:Projection:DataModel           -0.05      0.19    -0.42     0.32
## Motion:Shading:Projection:DataModel            0.27      0.19    -0.10     0.63
## Device:Motion:Shading:Projection:DataModel     0.20      0.37    -0.51     0.93
## shape_Intercept                                1.30      0.02     1.26     1.35
## shape_Device                                   0.22      0.05     0.13     0.31
## shape_Motion                                   0.32      0.05     0.22     0.41
## shape_Device:Motion                            0.20      0.10     0.01     0.40
##                                            Rhat Bulk_ESS Tail_ESS
## Intercept                                  1.00     1361     2099
## Device                                     1.00     3072     3488
## Motion                                     1.00     2404     3262
## Shading                                    1.00     3621     3624
## Projection                                 1.00     3976     3552
## DataModel                                  1.00     2930     3508
## Dimensionality                             1.00     3571     3331
## Device:Motion                              1.00     2564     3077
## Device:Shading                             1.00     3496     3657
## Motion:Shading                             1.00     3564     3428
## Device:Projection                          1.00     3704     3397
## Motion:Projection                          1.00     3707     3491
## Shading:Projection                         1.00     3915     3561
## Device:DataModel                           1.00     3009     3055
## Motion:DataModel                           1.00     2601     2912
## Shading:DataModel                          1.00     3854     3575
## Projection:DataModel                       1.00     3657     3362
## Device:Dimensionality                      1.00     3402     3159
## DataModel:Dimensionality                   1.00     3372     3664
## Device:Motion:Shading                      1.00     3559     3334
## Device:Motion:Projection                   1.00     3971     3807
## Device:Shading:Projection                  1.00     3918     3685
## Motion:Shading:Projection                  1.00     3751     3604
## Device:Motion:DataModel                    1.00     2701     3242
## Device:Shading:DataModel                   1.00     3919     3614
## Motion:Shading:DataModel                   1.00     3616     3661
## Device:Projection:DataModel                1.00     3600     3452
## Motion:Projection:DataModel                1.00     3756     3346
## Shading:Projection:DataModel               1.00     3820     3572
## Device:DataModel:Dimensionality            1.00     3543     3879
## Device:Motion:Shading:Projection           1.00     3723     3436
## Device:Motion:Shading:DataModel            1.00     3533     3542
## Device:Motion:Projection:DataModel         1.00     3781     3634
## Device:Shading:Projection:DataModel        1.00     3898     3767
## Motion:Shading:Projection:DataModel        1.00     3679     3492
## Device:Motion:Shading:Projection:DataModel 1.00     3725     3657
## shape_Intercept                            1.00     3592     3524
## shape_Device                               1.00     3810     3379
## shape_Motion                               1.00     3590     3237
## shape_Device:Motion                        1.00     3736     3499
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Now we have the models. Before we head towards guiding questions, we first a helper function for plotting the results.

draw_coefs <- function(m, 
                       target_coefs, 
                       title = "", 
                       fill_color = "white" , 
                       breaks_x = seq(0, 1, by = 0.5)
                       ){
  post_samples <- posterior_samples(m, pars = "b_") 
  
  b_CIs <-  post_samples %>% 
            select(target_coefs) %>% 
            gather(key = "coefs", value = "value") %>% 
            group_by(coefs) %>%
            median_qi() %>%
            ungroup() %>%
            mutate(lower95 = exp(.lower), upper95 = exp(.upper), median = exp(value))%>%
            select(coefs, median, lower95, upper95) %>%
            mutate(coefs = factor(coefs, levels = rev(target_coefs)), 
                 s_lower95 = round(lower95, 3), s_upper95 = round(upper95, 3),
                 CIs = paste0(round(median, 3), ' [', s_lower95, ',',  s_upper95, ']'))
    
  b_samples <-  post_samples %>% 
    select(target_coefs) %>% 
    gather(key = "coefs", value = "value") %>% 
    ungroup() 
  
 pr <- b_samples %>% 
    mutate(`Pr > 0` = if_else(value >= 0 , 1, 0)) %>%
    group_by(coefs) %>%
    summarise(`Pr > 0` = sum(`Pr > 0`) / nrow(post_samples)) %>%
    mutate(`Pr > 0` = round(`Pr > 0` * 100, 2))
    
 
  g <- b_samples %>%
    mutate(value = exp(value)) %>%
    ggplot(aes(x = value, y = coefs)) + 
      geom_richtext(data = b_CIs, mapping = aes(x = tail(breaks_x, n = 1) * 1.1, y = coefs, label = CIs), 
                 size = 2, 
                 fill = NA, label.color = NA, # remove background and outline
                 label.padding = grid::unit(rep(0, 4), "pt")) + 
      geom_richtext(data = pr, mapping = aes(x = tail(breaks_x, n = 1) * 1.25, y = coefs, label = `Pr > 0`), 
                 size = 2, 
                 fill = NA, label.color = NA, # remove background and outline
                 label.padding = grid::unit(rep(0, 4), "pt")) + 
      geom_density_ridges(
                  fill = fill_color, color = fill_color, 
                  rel_min_height = 0.0001,
                  alpha = 0.5,
                  scale = .9,
                  size = 0.5
                  ) + 
      geom_vline(xintercept = 1, linetype = "dashed", size = LINE_SIZE) +
      scale_x_continuous(breaks = breaks_x, expand = expansion(mult = c(0.01, .5))) +
      coord_cartesian(xlim = c(breaks_x[1], tail(breaks_x, n = 1))) +
      ggtitle(title) + xlab("exp(beta)") + ylab("") 
  
  return(g)
}

Results

GQ1: How do the six manipulated variables affect participants’ performance?

To answer the first guiding question, we look at the model coefficients of each experimental variable from fixed effects. The coefficients from fixed effects tell us the overall weight of each variable when conditioning on just that variable for an average; note that the sign and direction of an effect depends on coding.

target_coefs_rq1 <- c("b_Device", "b_Motion", "b_Projection", "b_Shading", "b_Dimensionality", "b_DataModel")


g1_rq1 <- draw_coefs(m_brms_error_hurdle_gamma, 
                 target_coefs_rq1, 
                 "Error magnitude", 
                 ERROR_COLOR, 
                 breaks_x = seq(0.75, 1.5, by = 0.25))
g2_rq1 <- draw_coefs(m_brms_time_gamma, 
                 target_coefs_rq1, 
                 "Response time",
                 TIME_COLOR, 
                 breaks_x = seq(0.85, 2.2, by = 0.45))

g_combined_rq1 <- ggarrange(g1_rq1, g2_rq1, ncol = 2)

ggsave(g_combined_rq1, filename = "pic/RQ1-mean-effects.pdf", width = 8, height = 2.25, units = "in")

g_combined_rq1

GQ2: How do the experimental variables interact with each other?

Similar to our investigation for the first guiding question, we use the model coefficients to investigate the interactions effects between variables. We only considered two-way interaction effects for simplicity; more complicated interaction effects were implied by the analyses below in GQ4.

target_coefs_rq2 <- c("b_Device:Motion", 
                  "b_Device:Shading", 
                  "b_Device:Projection",
                  "b_Device:Dimensionality", 
                  "b_Device:DataModel", 
                  "b_Motion:Shading", 
                  "b_Motion:Projection", 
                  "b_Motion:DataModel", 
                  "b_Shading:Projection", 
                  "b_Shading:DataModel", 
                  "b_Projection:DataModel", 
                  "b_DataModel:Dimensionality")


g1_rq2 <- draw_coefs(m_brms_error_hurdle_gamma, 
                 target_coefs_rq2, 
                 "Error magnitude", 
                 ERROR_COLOR, 
                 breaks_x = seq(0.5, 2.1, by = 0.5))
g2_rq2 <- draw_coefs(m_brms_time_gamma, 
                 target_coefs_rq2, 
                 "Response time",
                 TIME_COLOR, 
                 breaks_x = seq(0.5, 2.1, by = 0.5))

g_combined_rq2 <- ggarrange(g1_rq2, g2_rq2, ncol = 2)

ggsave(g_combined_rq2, filename = "pic/RQ2-interaction-effects.pdf", width = 8.5, height = 3.75, units = "in")

g_combined_rq2

GQ3: How does each manipulation affect participants’ performance?

To compare different levels of manipulations, we draw posterior samples from the model to get an estimate of mean for each measure. These posterior distributions show how an average participant does when conditioning on that level of manipulation.

target_order <- c("Device", "Motion", "Projection", "Shading", "Dimensionality", "DataModel")
m <- m_brms_error_hurdle_gamma
 draw_post <- function( m, 
                        title = "", 
                        fill_color = "white" , 
                        palette_in, 
                        adjust_in,
                        breaks_x = seq(0, 1, by = 0.5)
                        ){
   
  post_samples <- trials %>%
        # these use theoritical centers, not practical centers...
        # data_grid(Device = c(-0.5, 0.5), Motion = 0, Projection = 0, Shading = 0, Dimensionality = 0, DataModel = 0) %>%
        # add_row(trials %>%data_grid(Device = 0, Motion = c(-0.5, 0.5), Projection = 0, Shading = 0, Dimensionality = 0, DataModel = 0))  %>%
        # add_row(trials %>%data_grid(Device = 0, Motion = 0, Projection = c(-0.5, 0.5), Shading = 0, Dimensionality = 0, DataModel = 0))  %>%
        # add_row(trials %>%data_grid(Device = 0, Motion = 0, Projection = 0, Shading = c(-0.5, 0, 0.5), Dimensionality = 0, DataModel = 0))  %>%
        # add_row(trials %>%data_grid(Device = 0, Motion = 0, Projection = 0, Shading = 0, Dimensionality = c(-0.5, 0.5), DataModel = 0))  %>%
        # add_row(trials %>%data_grid(Device = 0, Motion = 0, Projection = 0, Shading = 0, Dimensionality = 0, DataModel = c(-0.5, 0.5)))  %>%
    select(Device, Motion, Projection, Shading, Dimensionality, DataModel) %>%
    distinct() %>%
        #View(post_samples)
    add_fitted_draws(m, value = "sample_value", scale = "response", re_formula = NA, allow_new_levels = F) %>%
     select(-c(".chain", ".iteration"))%>%
     ungroup()%>%
     select(Device, Motion, Projection, Shading, Dimensionality, DataModel, sample_value) %>%
     gather(key = "levels_names", value = "levels_values", Device, Motion, Projection, Shading, Dimensionality, DataModel)  %>%
     mutate(levels_names = factor(levels_names, levels = rev(target_order)), levels_values = factor(levels_values))

  b_CIs <-  post_samples %>%
            group_by(levels_names, levels_values) %>%
            median_qi()  %>%
            ungroup() %>%
            mutate(lower95 = .lower, upper95 = .upper, median = sample_value)%>%
            mutate(levels_values = factor(levels_values),
                 s_lower95 = round(lower95, 3), s_upper95 = round(upper95, 3),
                 CIs = paste0(round(median, 3), ' [', s_lower95, ',',  s_upper95, ']'))

  g <- post_samples %>%
      ggplot(aes(x = sample_value, y = levels_names)) +
      geom_density_ridges(
                aes(fill = NA, color = levels_values),
                rel_min_height = 0.0005,
                alpha = 0.5,
                scale = .8,
                size = 0.4
                ) +
      geom_richtext(data = b_CIs, mapping = aes(x = tail(breaks_x, n = 1) * 1.2,
                                                y = levels_names,
                                                color = levels_values,
                                                label = CIs),
                 label.color = NA,
                 size = 2,
                 fill = NA,
                 label.padding = grid::unit(rep(0, 4), "pt")) +
      #scale_x_continuous(breaks = breaks_x, expand = expansion(mult = c(0.01, .5))) +
      scale_fill_manual(values = palette_in) +
      scale_color_manual(values = palette_in) +
      #coord_cartesian(xlim = c(breaks_x[1], tail(breaks_x, n = 1))) +
      ggtitle(title) + xlab("") + ylab("")

   return(g)
 }

g1_rq3 <- draw_post(m_brms_error_hurdle_gamma, 
                 "Error magnitude (pt)", 
                 ERROR_COLOR, 
                 PALETTE_ERROR,
                 adjust_in = 1, 
                 breaks_x = seq(0, 10, by = 6))

g2_rq3 <- draw_post(m_brms_time_gamma, 
                 "Response time (s)",
                 TIME_COLOR, 
                 PALETTE_TIME,
                 adjust_in = 1, 
                 breaks_x = seq(0, 18, by = 12))

g_combined_rq3 <- ggarrange(g1_rq3, g2_rq3, ncol = 2)

ggsave(g_combined_rq3, filename = "pic/RQ3-div.pdf", width = 6, height = 3.4, units = "in")

g_combined_rq3

GQ4: Which resulting scatterplot designs are more effective for our task?

To investigate the effectiveness of each combination, we draw posterior samples from the model to get an estimate of all possible combinations. These posterior distributions show that how an average participant does with one combination, which corresponds to an experimental condition. Given our primary interest in error magnitude, we form a rank list of performance based on median of mean distribution.

 global_order <- NULL

 draw_all <- function( m, 
                        title = "", 
                        fill_color = "white" , 
                        palette_in, 
                        adjust_in,
                        breaks_x = seq(0, 1, by = 0.5),
                        which = "error"
                        ){
  
  post_samples <- trials %>%
    select(Device, Motion, Projection, Shading, Dimensionality, DataModel) %>%
    distinct() %>%
    add_fitted_draws(m, value = "sample_value", scale = "response", re_formula = NA, allow_new_levels = FALSE) %>%
    select(-c(".chain", ".iteration"))%>%
    ungroup()%>%
    select(Device, Motion, Projection, Shading, Dimensionality, DataModel, sample_value) %>%
    mutate(  # code it back so that we know which one is which...
             Device = recode(Device,  `0.5` = "HMD", `-0.5` = "Desktop"),
             Motion = recode(Motion,  `0.5` = "Dynamic", `-0.5` = "Static"),
             Projection = recode(Projection, `0.5` = "Persp.", `-0.5` = "Orth."),
             Shading = recode(Shading, `0` = "Simple", `0.5` = "A.O.", `-0.5` = "Flat"),
             DataModel = recode(DataModel, `0.5` = "Image Data", `-0.5` =  "Text Data"),
             Dimensionality = recode(Dimensionality, `-0.5` = "2D Embedding",  `0.5` = "3D Embedding"), 
             Condition = paste(Device, Motion, Projection, Shading, Dimensionality, DataModel, sep = seper))
  
  b_CIs <-  post_samples %>%
            select(Condition, sample_value) %>%
            group_by(Condition) %>%
            median_qi()  %>%
            ungroup() %>%
            mutate(lower95 = .lower, upper95 = .upper, median = sample_value, 
                   s_lower95 = round(lower95, 2), s_upper95 = round(upper95, 2),
                   CIs = paste0(round(median, 2), ' [', s_lower95, ',',  s_upper95, ']')) 
  if(which == "error"){
    b_CIs <- b_CIs %>%
            mutate(Condition = fct_reorder(Condition, .x = median, .desc = TRUE))
    global_order <<- levels(b_CIs$Condition)
  }else{
    b_CIs <- b_CIs %>%
            mutate(Condition = factor(Condition, levels = global_order))
  }
  
  
  median_max <- max(b_CIs$median)
  median_min <- min(b_CIs$median)
  
  g <- post_samples %>%
      ggplot(aes(x = sample_value, y = Condition)) +
      geom_vline(xintercept = median_max, linetype = "dashed", color = "black", size = LINE_SIZE) +
      geom_vline(xintercept = median_min, linetype = "dashed", color = "black", size = LINE_SIZE) +
      geom_richtext(data = b_CIs, mapping = aes(x = tail(breaks_x, n = 1) * 1.8, 
                                                y = Condition, 
                                                label = CIs), 
                 label.color = NA,
                 color = TEXT_COLOR,
                 size = 3, 
                 fill = NA, 
                 hjust = 1,
                 label.padding = grid::unit(rep(0, 4), "pt")) + 
      geom_density_ridges(
                  fill = fill_color,
                  color = fill_color,
                  rel_min_height = 0.00005,
                  alpha = 0.5,
                  scale = .9,
                  size = 0.5
                  ) +
      scale_x_continuous(breaks = breaks_x, expand = expansion(mult = c(0.01, .8))) +
      coord_cartesian(xlim = c(breaks_x[1], tail(breaks_x, n = 1))) +
      ggtitle(title) + xlab("") + ylab("") 
  g
  return(g)
}

g1_rq4 <- draw_all(m_brms_error_hurdle_gamma, 
                 "Error magnitude (pt)", 
                 ERROR_COLOR, 
                 PALETTE_ERROR,
                 adjust_in = 1, 
                 breaks_x = seq(0, 18, by = 6),
                 which = "error")

g2_rq4 <- draw_all(m_brms_time_gamma, 
                 "Response time (s)",
                 TIME_COLOR, 
                 PALETTE_TIME,
                 adjust_in = 1, 
                 breaks_x = seq(0, 36, by = 12),
                 which = "time")

g_combined_rq4 <- ggarrange(g1_rq4, g2_rq4, ncol = 2)

ggsave(g_combined_rq4, filename = "pic/RQ4-all.pdf", width = 14, height = 15, units = "in")
ggsave(g1_rq4, filename = "pic/RQ4-error.pdf", width = 8, height = 12.5, units = "in")
ggsave(g2_rq4, filename = "pic/RQ4-time.pdf", width = 8, height = 12.5, units = "in")

g_combined_rq4

This is the encoding to design the bubbles maps.

 draw_tiles <- function( m, 
                        title = "", 
                        fill_color = "white" , 
                        stop_color = "white",
                        palette_in, 
                        adjust_in,
                        breaks_x = seq(0, 1, by = 0.5),
                        which = "error"
                        ){
  
  post_samples <- trials %>%
    select(Device, Motion, Projection, Shading, Dimensionality, DataModel) %>%
    distinct() %>%
    add_fitted_draws(m, value = "sample_value", scale = "response", re_formula = NA, allow_new_levels = FALSE) %>%
    select(-c(".chain", ".iteration"))%>%
    ungroup()%>%
    select(Device, Motion, Projection, Shading, Dimensionality, DataModel, sample_value) %>%
    mutate(  Device = recode(Device,  `0.5` = "2 HMD", `-0.5` = "1 Desktop"),
             Motion = recode(Motion,  `0.5` = "2 Dynamic", `-0.5` = "1 Static"),
             Projection = recode(Projection, `-0.5` = "1 Orthographic", `0.5` = "2 Perspective"),
             Shading = recode(Shading, `-0.5` = "1 Flat Shading", `0` = "2 Simple", `0.5` = "3 Ambient Occlusion"),
             DataModel = recode(DataModel, `0.5` = "2 Image Data", `-0.5` =  "1 Text Data"),
             Dimensionality = recode(Dimensionality, `-0.5` = "2D Embedding",  `0.5` = "3D Embedding"), 
             Condition = paste(Dimensionality, Device, Motion, DataModel, Projection, Shading,  sep = seper))
  
  b_CIs <-  post_samples %>%
            select(Condition, sample_value) %>%
            group_by(Condition) %>%
            median_qi()  %>%
            ungroup() %>%
            mutate(lower95 = .lower, upper95 = .upper, median = sample_value, 
                   s_lower95 = round(lower95, 3), s_upper95 = round(upper95, 3),
                   CIs = paste0(round(median, 3), ' [', s_lower95, ',',  s_upper95, ']'),
                   index = row_number()
                   ) 
  scale_factor <- 2.8
  g <- b_CIs %>%
      ggplot() +
      theme(legend.position = "right", panel.grid.major = element_blank()) +
      coord_fixed(ratio = 1) +
      geom_tile(aes(x = 47, y = scale_factor * index), width = scale_factor, height = scale_factor, size = 0.1, fill = NA, color = "black") + 
      geom_circle(aes(color = median, x0 = 47, y0 = scale_factor * index, r = adjust_in * sqrt(upper95)), fill = NA, size = 0.1) +
      geom_circle(aes(fill = median, x0 = 47, y0 = scale_factor * index, r = adjust_in * sqrt(median)), color = NA,  size = 0.1) +
      geom_circle(aes(x0 = 47, y0 = scale_factor * index, r = adjust_in * sqrt(lower95)), color = "white", fill = NA, size = 0.1) +
      geom_text(aes(x = 20, y = scale_factor * index, label = Condition), size = 2) +
      scale_x_continuous(limits = c(0, 50)) +
      scale_fill_gradient(low = stop_color, high = fill_color, limits = c(breaks_x[1], tail(breaks_x, n = 1)), guide = "colourbar") +
      scale_color_gradient(low = stop_color, high = fill_color, limits = c(breaks_x[1], tail(breaks_x, n = 1)), guide = "colourbar") + 
      scale_y_reverse() + 
      ggtitle(title) + xlab("") + ylab("") 
  g
  return(g)
}

g1_rq4_tiles <- draw_tiles(m_brms_error_hurdle_gamma, 
                 "Error magnitude (pt)", 
                 ERROR_COLOR, 
                  "white",
                 PALETTE_ERROR,
                 adjust_in = 0.35, 
                 breaks_x = seq(3, 15, by = 4),
                 which = "error")

g2_rq4_tiles <- draw_tiles(m_brms_time_gamma, 
                 "Response time (s)",
                "#699923", #8bba45 
                "#f0fcd7", 
                 PALETTE_TIME,
                 adjust_in = 0.24, 
                 breaks_x = seq(0, 36, by = 12),
                 which = "time")

g_combined_rq4_tiles <- ggarrange(g1_rq4_tiles, g2_rq4_tiles, ncol = 2)

ggsave(g_combined_rq4_tiles, filename = "pic/RQ4-bubbles.pdf", width = 10.5, height = 19.5, units = "in")

g_combined_rq4_tiles

Misc figures

Demographics

demographics <- metas %>% select(familiarML, familiarDR, familiarVIS, familiarVR)
demographics_stack <- stack(demographics)
demographics_stack$overall <- paste(demographics_stack$values, demographics_stack$ind, sep = " ")

demographics_table <- aggregate(data = demographics_stack, values ~ overall, length)
demographics_table$item <- sapply(demographics_table$overall, function(x){return(strsplit(x, split = " ")[[1]][2])})  
demographics_table$option <- sapply(demographics_table$overall, function(x){return(strsplit(x, split = " ")[[1]][1])})  

df <- data.frame(
  item = c("normal", "glasses", "contacts", "others", "a little bit", "not at all"),
  values = c(10, 5, 5, 2, 9, 23)
)

demo_color <- "#6959a6"

g_demographics <- ggplot() +
  theme(legend.position = "top", panel.grid.major = element_blank()) +
  #geom_tile(data = demographics_table, aes(x = option, y = item, fill = values)) +
  geom_point(data = demographics_table, aes(color = values, x = as.numeric(option), y = item, size = sqrt(values)), stroke = 1) +
  geom_point(data = df, aes(color = values, x = 8, y = item, size = sqrt(values)), stroke = 1) +
  scale_fill_gradient(low = "#ffffff", high = demo_color, limits = c(0, 25), guide = "colourbar") +
  scale_color_gradient(low = "#ffffff", high = demo_color, limits = c(0, 25), guide = "colourbar") + 
  scale_x_continuous(limit = c(0, 9), breaks = seq(1, 7, by = 1)) + 
  xlab("scale") + ylab("")

ggsave(g_demographics, filename = "pic/g_demographics.pdf", width = 3.5, height = 3.4, units = "in")

g_demographics

Accuracy range examples

acc1 <- c(0.3978, 0.4054, 0.5998, 0.602, 0.605, 0.6052, 0.6528, 0.653, 0.6614, 0.6762, 0.6956, 0.7044, 0.713, 0.7162, 0.725, 0.73, 0.7326, 0.7362, 0.7398, 0.745, 0.7458, 0.7476, 0.751, 0.7538, 0.7668, 0.767, 0.7686, 0.7686, 0.769, 0.7744, 0.7748, 0.7754, 0.7758, 0.7778, 0.7788, 0.7794, 0.7806, 0.7828, 0.7834, 0.7848, 0.7858, 0.7882, 0.7888, 0.789, 0.7894, 0.7914, 0.7932, 0.7946, 0.7948, 0.7952, 0.7954, 0.7966, 0.7968, 0.7976, 0.799, 0.7992, 0.801, 0.8014, 0.8018, 0.8026, 0.8028, 0.8032, 0.804, 0.804, 0.8046, 0.805, 0.8056, 0.8058, 0.806, 0.8066, 0.8068, 0.8072, 0.8076, 0.808, 0.8084, 0.8098, 0.8108, 0.8108, 0.8108, 0.8118, 0.812, 0.8132, 0.8138, 0.8138, 0.8148, 0.816, 0.816, 0.8172, 0.8176, 0.8182, 0.8192, 0.8202, 0.8204, 0.8204, 0.8206, 0.8206, 0.8206, 0.8212, 0.8212, 0.822, 0.8224, 0.824, 0.8244, 0.8254, 0.8258, 0.8268, 0.8274, 0.8276, 0.8286, 0.8286, 0.829, 0.8302, 0.8314, 0.8314, 0.8326, 0.833, 0.8332, 0.8338, 0.834, 0.8352, 0.8358, 0.8364, 0.837, 0.8384, 0.8386, 0.8386, 0.8396, 0.8406, 0.8408, 0.8408, 0.8414, 0.842, 0.8428, 0.8436, 0.8442, 0.8444, 0.8444, 0.8446, 0.845, 0.8452, 0.847, 0.8472, 0.8472, 0.8474, 0.8474, 0.848, 0.85, 0.8504, 0.8506, 0.8514, 0.8516, 0.8518, 0.8534, 0.855, 0.8578, 0.858, 0.8582, 0.8584, 0.8584, 0.861, 0.8634, 0.8662, 0.8912, 0.8914, 0.8922, 0.8942, 0.8954, 0.896, 0.896, 0.896, 0.8962, 0.8964, 0.8966, 0.8968, 0.8972, 0.8974, 0.8978, 0.8978, 0.8978, 0.8982, 0.8982, 0.8984, 0.8984, 0.8986, 0.899, 0.8992, 0.8992, 0.8994, 0.8996, 0.8998, 0.8998, 0.8998, 0.8998, 0.8998, 0.9, 0.9, 0.9002, 0.9002, 0.9002, 0.9002, 0.9002, 0.9004, 0.9006, 0.9006, 0.9008, 0.901, 0.901, 0.9012, 0.9012, 0.9014, 0.9014, 0.9014, 0.9016, 0.9018, 0.9018, 0.9018, 0.902, 0.902, 0.902, 0.9022, 0.9024, 0.9026, 0.9026, 0.9026, 0.9026, 0.9026, 0.9028, 0.9028, 0.903, 0.903, 0.903, 0.903, 0.903, 0.9032, 0.9032, 0.9036, 0.9036, 0.9036, 0.9036, 0.9036, 0.9038, 0.9038, 0.9038, 0.9038, 0.9044, 0.9044, 0.9044, 0.9044, 0.9046, 0.9046, 0.905, 0.9052, 0.9052, 0.9052, 0.9054, 0.9054, 0.9056, 0.9056, 0.9056, 0.9058, 0.9058, 0.9058, 0.9058, 0.906, 0.906, 0.906, 0.906, 0.906, 0.9062, 0.9064, 0.9064, 0.9066, 0.9066, 0.9068, 0.9068, 0.9068, 0.9068, 0.907, 0.907, 0.9074, 0.9074, 0.9074, 0.908, 0.908, 0.9086, 0.9092, 0.9092, 0.9092, 0.9092, 0.9094, 0.9094, 0.9098, 0.9098, 0.91, 0.9106, 0.9114, 0.9114, 0.9126, 0.913, 0.9146)


acc2 <- c(0.192, 0.246, 0.394, 0.405, 0.427, 0.439, 0.476, 0.485, 0.488, 0.488, 0.49, 0.491, 0.491, 0.492, 0.492, 0.494, 0.496, 0.496, 0.497, 0.497, 0.497, 0.498, 0.498, 0.499, 0.499, 0.499, 0.5, 0.5, 0.5, 0.5, 0.5, 0.501, 0.504, 0.505, 0.505, 0.506, 0.506, 0.507, 0.507, 0.508, 0.508, 0.508, 0.508, 0.509, 0.509, 0.51, 0.51, 0.511, 0.511, 0.512, 0.512, 0.514, 0.515, 0.515, 0.515, 0.516, 0.522, 0.526, 0.53, 0.531, 0.531, 0.531, 0.534, 0.542, 0.548, 0.559, 0.564, 0.6, 0.625, 0.658, 0.681, 0.697, 0.722, 0.724, 0.726, 0.731, 0.735, 0.739, 0.739, 0.739, 0.739, 0.743, 0.747, 0.748, 0.753, 0.762, 0.764, 0.767, 0.776, 0.79, 0.79, 0.792, 0.795, 0.802, 0.803, 0.81, 0.813, 0.813, 0.817, 0.819, 0.82, 0.82, 0.823, 0.824, 0.824, 0.826, 0.828, 0.828, 0.834, 0.837, 0.838, 0.839, 0.841, 0.841, 0.842, 0.842, 0.843, 0.843, 0.844, 0.847, 0.848, 0.852, 0.853, 0.854, 0.854, 0.854, 0.856, 0.857, 0.857, 0.857, 0.859, 0.86, 0.862, 0.863, 0.863, 0.864, 0.865, 0.866, 0.871, 0.871, 0.874, 0.874, 0.875, 0.876, 0.879, 0.88, 0.881, 0.883, 0.883, 0.883, 0.884, 0.886, 0.886, 0.886, 0.887, 0.887, 0.888, 0.891, 0.891, 0.893, 0.893, 0.894, 0.896, 0.896, 0.896, 0.897, 0.898, 0.899, 0.901, 0.901, 0.901, 0.901, 0.902, 0.903, 0.904, 0.904, 0.904, 0.905, 0.905, 0.905, 0.907, 0.908, 0.908, 0.909, 0.91, 0.91, 0.91, 0.911, 0.911, 0.911, 0.911, 0.912, 0.912, 0.913, 0.913, 0.914, 0.917, 0.917, 0.92, 0.922)

acc3 <- c(0.838, 0.839, 0.8405, 0.8505, 0.8575, 0.8575, 0.865, 0.872, 0.872, 0.873, 0.8735, 0.8745, 0.882, 0.882, 0.8835, 0.885, 0.8865, 0.887, 0.887, 0.892, 0.892, 0.8935, 0.895, 0.8955, 0.896, 0.897, 0.8975, 0.898, 0.898, 0.8985, 0.8985, 0.899, 0.8995, 0.8995, 0.8995, 0.9, 0.9015, 0.902, 0.9025, 0.9035, 0.904, 0.9045, 0.9045, 0.905, 0.9055, 0.9055, 0.906, 0.906, 0.9075, 0.908, 0.9085, 0.9085, 0.9085, 0.9085, 0.9095, 0.9095, 0.9105, 0.9105, 0.9105, 0.911, 0.9115, 0.9115, 0.9115, 0.9115, 0.9125, 0.9125, 0.9125, 0.9125, 0.9125, 0.9125, 0.913, 0.9135, 0.9135, 0.9135, 0.914, 0.914, 0.914, 0.9145, 0.915, 0.915, 0.915, 0.9155, 0.9155, 0.9155, 0.9155, 0.916, 0.916, 0.916, 0.9165, 0.9165, 0.917, 0.917, 0.917, 0.917, 0.9175, 0.918, 0.918, 0.918, 0.918, 0.918, 0.918, 0.918, 0.9185, 0.9185, 0.9185, 0.9185, 0.9185, 0.919, 0.919, 0.919, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.92, 0.92, 0.92, 0.92, 0.9205, 0.9205, 0.9205, 0.9205, 0.921, 0.921, 0.9215, 0.9215, 0.9215, 0.9215, 0.9215, 0.922, 0.922, 0.922, 0.922, 0.922, 0.922, 0.9225, 0.9225, 0.9225, 0.9225, 0.923, 0.923, 0.923, 0.923, 0.9235, 0.9235, 0.9235, 0.9235, 0.9235, 0.9235, 0.924, 0.924, 0.924, 0.924, 0.9245, 0.9245, 0.9245, 0.9245, 0.9245, 0.9245, 0.9245, 0.9245, 0.925, 0.925, 0.925, 0.925, 0.9255, 0.9255, 0.9255, 0.9255, 0.9255, 0.9255, 0.9255, 0.9255, 0.926, 0.926, 0.926, 0.926, 0.926, 0.9265, 0.9265, 0.927, 0.927, 0.927, 0.927, 0.927, 0.927, 0.927, 0.927, 0.927, 0.9275, 0.9275, 0.9275, 0.9275, 0.928, 0.928, 0.928, 0.928, 0.928, 0.928, 0.928, 0.928, 0.9285, 0.9285, 0.9285, 0.9285, 0.9285, 0.9285, 0.929, 0.929, 0.929, 0.929, 0.929, 0.929, 0.9295, 0.9295, 0.93, 0.93, 0.93, 0.93, 0.93, 0.9305, 0.9305, 0.9305, 0.931, 0.9315, 0.9315, 0.9315, 0.9315, 0.932, 0.9325, 0.9325, 0.933, 0.9335, 0.9335, 0.9335, 0.934, 0.9345, 0.9345, 0.936, 0.936, 0.937, 0.937, 0.94, 0.9405)

bars_color <- "#6959a6"

g_bar1 <- ggplot() +
           theme(panel.grid.major.x = element_blank()) +
           geom_linerange(aes(ymax = 1, ymin = 0, x = acc1), size = 0.07, color = bars_color) +
           scale_x_continuous(limits=c(0,1), breaks = seq(0, 1, by = 0.1), expand=c(0,0)) +
           scale_y_continuous(limits=c(0,1), breaks = seq(0, 1, by = 1), expand=c(0,0)) +
           ggtitle("CIFAR-10")

g_bar2 <- ggplot() +
            theme(panel.grid.major.x = element_blank()) +
            geom_linerange(aes(ymax = 1, ymin = 0, x = acc2), size = 0.07, color = bars_color) +
            scale_x_continuous(limits=c(0,1), breaks = seq(0, 1, by = 0.1), expand=c(0,0)) +
            scale_y_continuous(limits=c(0,1), breaks = seq(0, 1, by = 1), expand=c(0,0))+
            ggtitle("bABI")

g_bar3 <- ggplot() +
            theme(panel.grid.major.x = element_blank()) +
            geom_linerange(aes(ymax = 1, ymin = 0, x = acc3), size = 0.07, color = bars_color) +
            scale_x_continuous(limits=c(0,1), breaks = seq(0, 1, by = 0.1), expand=c(0,0)) +
            scale_y_continuous(limits=c(0,1), breaks = seq(0, 1, by = 1), expand=c(0,0)) +
            ggtitle("F-MNIST")

g_bars <- ggarrange(g_bar1, g_bar2, g_bar3, ncol = 1, align="v")

ggsave(file = "pic/acc_bar.pdf", g_bars , width = 4, height = 2, units = "in")

g_bars